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ABSTRACT 


A  repairement  assumption  called  instantaneous  resurrection 
is  introduced  and  it  is  shown  that,  under  this  assumption,  the 
number  of  failures  of  an  equipment  in  time  t  is  a  Poisson  process 
with  leading  function  -log  R(t),  where  R(t)  is  the  reliability 
function  of  the  equipment.  This  fact  Is  used  in  the  estimation 
of  parameters  of  prior  distributions  via  the  appropriate  Poisson 
mixtures.  The  case  of  the  two-way  truncated  Gamma  prior  and  its 
limiting  forms  is  handled  in  some  detail.  A  preliminary  form  of 
a  goodness  of  fit  test  in  the  case  of  varying  operational  time 
in  attribute  failure  data  is  presented. 
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EVALUATION 


The  work  described  in  the  report  extends  previous  work  by 
Dr.  Bolis  which  provided  guides  to  structuring  Bayesian  Reliability 
Test  Plans  assuming  an  inverted  gamma  distribution  of  the  prior 
information  and  a  constant  failure  rate  of  the  tested  equipment.  The 
report  considers  cases  where  these  assumptions  do  not  hold.  The 
previous  work  is  currently  being  incorporated  into  a  standard  for 
Bayesian  Reliability  Tests.  The  results  presented  in  this  report  will 
be  used  for  modifying  the  standard. 


ANTHONY  COPPOLA 
Project  Engineer 


0.  INTRODUCTION 


In  this  report  we  introduce  a  repairment  assumption,  which  we  call 
resurrection.  Under  this  assumption  when  an  equipment  fails  and  is 
repaired,  the  repair  does  not  bring  the  equipment  to  a  brand  new  stage 
but  to  a  stage  of  a  functioning  equipment  of  the  same  type  and  age  that 
had  no  failures.  This  assumption  may  be  deemed  as  a  simplication  of  the 
actual  circumstances,  but  to  our  opinion,  it  is  a  lot  more  preferable 
than  the  renewal  assumption. 

We  show  in  theorem  1.2.1  that,  under  the  instantaneous  resurrection 
assumption,  the  number  of  failures  of  an  equipment  in  time  t  is  a  non- 
homogeneous  Poisson  process  with  leading  function  -  log  R(t),  where  R(t) 
is  the  reliability  function  of  the  equipment.  If  we  assume  a  suitable 
form  of  a  prior  distribution  of  the  parameters  of  the  reliability  function, 
the  resulting  Poisson  mixture  (marginal  distribution)  allows  the  estimation 
of  the  prior  distribution  from  prior  data.  When  the  prior  is  determined, 
Bayesian  acceptance  sampling  plans  can  be  constructed.  The  above  form 
the  content  of  Section  1 . 

In  Section  2  we  consider  the  case  of  constant  failure  rate,  where 
the  resurrection  and  renewal  assumptions  coincide.  The  two-way  truncated 
inverted  Gamma  distribution  as  well  as  all  of  its  limiting  forms  (truncated 
Pareto,  uniform,  etc.)  is  assumed  and  the  methods  of  estimating  its 
parameters  (mainly  the  method  of  moments)  are  presented. 

In  Section  3  we  indicate  how  the  methods  of  Section  2  can  be  carried 
out  (with  minor  adjustments)  to  the  case  of  variable  failure  rate.  As 
an  example,  the  Weibuil  time  to  failure  is  considered. 

In  Section  4  we  construct  a  statistic* akin  to  the  statistic,  to 
test  goodness  of  fit  in  the  case  the  failure  data  do  not  exhibit  the 
same  total  time  of  operation.  The  exact  mean  and  variance  of  this  statis¬ 
tic  as  well  as  the  cumulants  of  its  asymptotic  distribution  are  computed. 

We  do  not  know  how  these  quantities  are  affected  when  parameters  are 
estimated  from  thesample used  in  the  construction  of  this  statistic. 

This  report  is  of  a  theoretical  nature.  It  is  of  interest  to  check 
the  resurrection  assumption  by  using  failure  data  of  repairable  equipment 
with  known  distribution  of  time  to  failure  other  than  exponential.  Also, 
it  would  be  of  interest  to  see  the  amount  of  Improvement  of  fit  if  the 
two-way  truncated  inverted  Gamma  prior  distribution  is  used  on  the  mean 
time  to  failure  of  electronic  equipment  in  comparison  with  the  use  of 
the  full  inverted  Gamma.  Further,  the  question  of  settling  the  behaviour 
of  the  statistic  S  of  Section  4  when  parameters  are  estimated  is  important. 
Some  future  effort  may  be  dedicated  to  these  questions. 


1 .  GENERAL  THEORY 


1.1  THE  DISTRIBUTION  OF  THE  TIME  TC  FAILURE.  Let  f(t),  t>o  be  the 
density  function  of  the  time  to  first  failure  of  an  equipment.  We 
assume  that  the  function  f  is  continuous.  The  function 


R(t)  =  /  f(x)dx,  t>o, 

is  called  the  reliability  function  of  the  equipment  under  discussion. 
It  is  the  probability  that  the  equipment  will  first  fail  after  time  t 
The  failure  rate  X(t)  (or  hazard  rate,  or  force  of  mortal ity)  of  this 
equipment  is  defined  by  — -  — ~ 


X(t)  *  — -7-  log  R(t)  =  f(t)/R(t). 


Another  measure  of  reliability  is  the  mean  time  to  failure.  In  the  general 
case  of  an  aging  repairable  equipment,  it  is  more  precise  to  call  this 
measure  the  mean  time  to  first  failure.  If  lim  tR(t)  =  0,  t-*-«,  this 
mean  is  given  by 


(1 .1 .1 )  0  s  /°°  t  f(t)dt  =  /  R(t)dt. 

1  0  0 

1.2  THE  DISTRIBUTION  OF  ATTRIBUTE  FAILURE  DATA.  When  a  repairable  equipment 
is  put  to  life  testing  and  a  failure  occurs,  two  simple  and  opposite 
assumptions  can  be  made  with  regard  to  its  state  after  repair:  The  repair 
either  brought  the  equipment  to  the  start  of  a  brand  new  one  (renewal 
assumption)  or  to  the  start  of  an  equipment  of  a  failure- free  equipment  of 
the  same  age  (resurrection  assumption).  The  first  case  has  been  extensively 
studied  in  the  literature.  Although  the  actual  situation  lies  somewhere 
between  these  two  extreme  cases,  we  believe  that  the  resurrection  assumption 
is  more  realistic  in  the  case  of  nonconstant  failure  rate. 

When  a  repairable  equipment  is  tested,  we  shall  assume  instantaneous 
repair,  i.e.,  either  the  repair  time  is  negligible  or  it  is  not  included 
as  test  time.  Two  quantities  anter  in  such  testing:  the  test  time  T 
and  the  number  of  failures  r  during  the  time  T.  We  can  either  fix  the 
test  time  T  in  advance  and  observe  the  random  variable  r  (type  I  censoring) 
or  test  the  equipment  until  a  specified  number  of  failures  occur  (type  II~ 
censoring) .  While  a  type  I  censored  test  has  obvious  advantages  (test 
time  and  cost  well  controlled),  a  type  II  censoring  has  theoretical  as  well 
as  practical  interest.  For  example,  if  only  r  repairs  are  possible,  then 
one  would  be  interested  in  the  distribution  of  the  waiting  time  for  the 
(r+1)  th  failure. 

1.2.1  TYPE  I  CENSORING  UNDER  THE  INSTANTANEOUS  RESURRECTION 
ASSUMPTION.  It  is  surprising  that  the  distribution  of  the  number  of 
failures  r  in  constant  time  T  under  the  instantaneous  resurrection 
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assumption  is  always  a  Poisson  distribution  with  mean  -  log  R(T).  We 
formalize  and  prove  this  result  in  the  following: 

THEOREM  1.2.1.  Under  the  instantaneous  resurrection  assumption,  the 
di  stri buti on  of  the  number  of  fail  ures  of  ari  equipment  wi th  rel  iabJM  ity 
function  R(  t)  in  a_  f i xed  time  T  j_s_  the  Poi  sson  di  stri bution  wi  th  mean 

-  log  R(T) . 

PROOF:  The  density  function  of  the  time  tj  of  the  first  failure  is  of 
course  f(t-j),  tj>o.  The  cumulative  distribution  function  of  the  time  t2  of 
the  second  failure  is  given  by 

Prob  (t^x | t2>t-| }  =  Prob  {t-|<t2;ix}  /  Prob  {t2>t-|}  2  (F(x)  -  F(t-| )  )/R(t^) , 
x>t-| , 

where  F  denotes  the  cumulative  distribution  function  of  the  time  to 
failure.  At  this  point  we  used  the  instantaneous  resurrection  assumption. 
Thus,  the  density  function  of  t2  is  f (t2)/R(t-j ) ,  t2>t-| .  Similarly,  the 
density  function  of  the  time  t j  of  the  ith  failure  is  f (tj )/R(t1-_] ) , 
t-pt-j^i,  l<i<r.  Finally,  the  probability  that  no  failures  occurred  between 
times  tf  and  T  is  R(T)/R(tr).  Thus  the  joint  density  function  of  the  times 
t-j ,  t2,  , ...  tr  with  0<tT<t2<-..<triT  is 


f(t-j)  f(t2) 
R(  t-j ) 


lUd  Mil 

R(tr_i )  R(tr) 


R(T)  ,\  (t-j)  ...  X  (tr) 


Therefore,  the  probability  that  r  failures  occurred  in  fixed  time  T  is 


T 

T 

T 

PT(r)  =  R(T) 

f 

/ 

...  /  X(ti 

tr-l 

0 

*1 

T 

T 

T 

=  RO) 

r 

J 

f 

...  /  X(t,) 

0 

tl 

tr-2 

R(T) 

logr 

( 1  / R  ( T )),  r-o,  1 

r! 


)  X(t2 )  • . .  X(tr)  dtr  ...  dt-j 

—  X(tr  l)  log  ^Vl j,  dtr  i  ...dt 
R(T) 

,  » •  *  q.e.d. 


More  generally,  let  E(t)  be  the  number  of  failures  of  the  equipment  in 


time  t.  By  reasoning  as  in  theorem 

(1  .2.1 )  Prob  (E(t)  -  £(s)  =  r} 

i.e.  the  stochastic  process  ^(t)  is 
leading  function  -log  R ( t )  (cf.  page 


1,2,1,  we  obtain 


=  pr  1ogr  (R(s)/R(t))  ,t>s, 

a  non-homogeneous  Poisson  process  with 

46  of  (1  ;) 
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1.2.2.  TYPE  II  CENSORING  UNDER  INSTANTANEOUS  RESURRECTION.  Let  T 
denote  the  waiting  time  until  the  rth  failure  (r  fixed).  We  shall  show 
that  the  random  variable  -  log  R(T)  has  the  Gamma  distribution  with  shape 
parameter  r  and  scale  parameter  1.  In  view  of  theorem  1.2.1  this  does 
not  seem  surprising. 

THEOREM  1.2.2.  Under  the  instantaneous  resurrection  assumption ,  the 

density  function  of  the  waiting  time  T_  unti  1  the  rth  fai  1  ure  is  gi_ven 
~  "  '  “  . . 

(1.2.2)  fr(T)  s  |“-|  j ,  f(T)  log^1  (1/R(T) ) ,  T>o. 

PROOF:  By  reasoning  along  lines  similar  to  those  of  the  proof  of 
theorem  1.2.1  we  obtain  that  the  cumulative  distribution  function  of  T 
is  given  oy 


X  X 

Prob  {T<x}  =  /*  /  •■•  /  X(ti)  ...  \(t  , )  f(t  )  dt  ...  dt,  . 

—  oo  tr_i  '  r-i  r  r  i 


[■'j  differentiating  this  function  with  respect  to  x  and  replacing  x  by  T 
we  obtain  (1.2.2).  q.e.d. 


The  average  waiting  time  to  the  rth  failure  will  be  called  the  mean  time 
to  rth  failure  and  it  will  be  denoted  by  Or. 

Thus , 

(1.2.3)  0  5  /*  T  logr_1  (1/R(T))  f(T)  dT  r*l ,  2,  ... 

(r-1 ) !  0 


This  agrees  with  the  formula  (1.1.1). 


1.2.3  THE  INSTANTANEOUS  RENEWAL  ASSUMPTION.  The  distribution  of 
attribute  failure  data  under  the  instantaneous  renewal  assumption  is 
well  worked  out  in  the  literature  (cf.  e.g.  (1),  (2),  (3)).  In  the 
type  I  situation  we  have 

(1.2.4)  PT(r)  =  /T  /T  ...  /T  f(t.)  )-R(T-t_ )  dt  ...  dt, 

T  0  t1  tr.!  12  1  r  r  1 


whereas  in  the  type  II  case  we  obtain 

J 


(1  .2.5)  f  (T)  =  fJ  fT 

'  r\ 


'1 


/;  f(t1)  f(t2-t1)  •••  f(T-tr.1} 


dt 


r-1 


dt. 


T>o ,  r=2,  3,  ...  , 


f](T)  *  f  (T) . 
If  we  set  F0(T)  =  1,  Fp(T) 


fr (t)dt,  r  =  1 , 


2, 


(T>o ) ,  we  get 


4 


PT(r)  -  Fr(T)  -  Fr+,  (T). 


In  general,  the  integrations  occurring  in  (1.2.4)  and  (1.2.5)  are  difficult 
to  carry  out  in  closed  form.  It  is  well  known  that  the  instantaneous 
resurrection  assumption  Is  equivalent  to  the  instantaneous  renewal  assumption 
if  and  only  If  the  distribution  of  the  time  to  failure  Is  exponential. 

1.2.4  CENSORING  WITHOUT  REPLACEMENT.  Suppose  that  n  identical  equip¬ 
ments  are  put  to  life  testing  until  a  beforehand  specified  time  T.  Suppose 
further  that  r  equipments  failed  at  times  tj  t2  ^  . . .  <_  t  T.  The 
joint  density  of  these  times  Is  given  hy  r 

Hi  f(t-, )  ...  f(tr)  R(T)n‘r  o<  t,  <  ...  <tr<T 
(n-r)! 

and  the  probability  function  of  the  numer  of  failures  r  is  binomial: 

(1.2.6)  PT(r)  *  (J)  (1 -R(T) )r  R(T)n_r  r  -  0,  1 . n. 

Similarly, if  n  equipments  are  tested  until  the  rth  failure  occurs  (when  an  equip- 
fails.it  is  discarded) .the  times  of  failure  0  ^  tr  j  <  t  s  T  have  joint 

density  function 

Hi  fft,)  ...  f(t  -,)  f (T)  R(T)n'r,  o  <  t1  <  ...  <  tr_-|  <  1 
(n-r)! 


anu  the  density  function  of  T  is  given  by 

(1.2.7)  fr(T)  -  r  (J)  f  (T)  (l-R(T))r'1  R(T)n'r.  T>o 

i.e.  the  random  variable  R(T)  has  the  Beta  distribution  with  parameters 
n-r+1  and  r. 

1.3  BAYESIAN  ESTIMATION  TECHNIQUES  OF  RELIABILITY  CHARACTERISTICS.  Usually 
the  time  to  failure  distribution  depends  on  a  (possibly  vector-valued) 
parameter  a,  i.e.  fit)  is  of  the  form 

f(t|a),  t>0,  acG] ^ 

where  Gi  is  the  range  of  the  parameter  a  (usually  a  nice  subset  of  the 
Euclidean  space  R^).  Knowledge  of  the  parameter  a  entails  complete 
knowledge  of  the  reliability  characteri sti cs  of  the  equipment.  Classical 
estimation  techniques  apply  to  failure  data  of  the  different  kinds  mentioned 
In  the  previous  section  1.2  to  yield  point  estimates  and  confidence  regions 
for  ot.  Also  tests  of  hypothesis  (reliability  demonstration)  can  be  devised 
for  a.  These  latter  tests  can  be  acceptance  tests  with  their  power  functions 


delineating  the  producer's  and  consumer's  risks.  Here,  we  confine  our¬ 
selves  to  Bayesian  techniques. 

1.3.1  DETERMINATION  OF  THE  PRIOR  DISTRIBUTION.  The  choice  of  the 
form  of  the  prior  distribution  of  a  is  made  on  the  basis  of  computational 
tractability,  sufficient  flexibility  to  accommodate  a  wide  variety  of 
distribution  shapes  and,  if  possible,  on  whether  it  gives  rise  to  conjugate 
pairs,  i.e.,  on  whether  the  posterior  distribution  of  a  given  the  failure 
data  will  be  of  the  same  form  as  the  prior.  Once  the  form  of  the  prior  has 
been  decided  upon,  it  remains  to  determine  its  parameters.  We  shall  assume 
that  prior  failure  data  of  the  equipment  in  question  are  available.  Since 
a  is  unobservable,  the  parameters  of  the  prior  can  be  estimated  from  the 
marginal  distribution  of  the  available  failure  data,  provided  that  this 
marginal  distribution  is  uniquely  determined  by  the  prior.  In  this  case, 
we  say  that  the  prior  distribution  is  identi fiabl e  by  the  marginal  distri¬ 
bution  of  the  data. 

Let  us  be  more  specific.  Let 
g(a;8),  aeG]CR  8cG2«R m 

be  the  prior  density  of  a,  where  8  is  a  (possibly  vector-valued)  parameter 
to  be  determined  (estimated).  If  the  available  failure  date  are  attribute 
data  of  the  type  I  censoring,  then  the  marginal  distribution  in  question 
will  be 

(1.3.1)  /  Py  ( r | a)  g(a;S)  da 

G1 

where  Py(r|a)  is  given  either  by  theorem  1.2.1  or  by  (1.2.4)  or  by  (1.2.6) 
depending  on  the  nature  of  the  available  data.  The  parameter  8  can  be 
estimated  by  classical  methods  under  the  hypothesis  of  identifi ability. 

The  appropriateness  of  the  prior  can  be  checked  by  a  x.2  -  goodness  of  fit 
test.  If  the  data  are  of  the  type  II  censoring  (a  rare  situation),  a 
similar  procedure  can  be  applied. 


The  procedure  described  so  far  in  this  section  was  initiated  by 
Schafer  et.al.  (4)  for  the  case  of  exponential  time  to  failure  with  an 
inverted  Gamma  prior  distribution  for  the  mean  time  to  failure.  The 
method  of  moments  was  used  In  (4)  whereas  Goel  (5)  used  the  maximum  likeli¬ 
hood  method  in  the  same  situation.  When  the  available  failure  data  are 


of  the  form  (r^,  T. ),  i*l,  2,  ...  n  (the  ith  equipment  had  r-j  failures  in 
time  Tj)  it  was  proposed  by  us  in  (6)  that  in  order  to  estimate  the 
parameter(s)  of  the  prior  distribution,  the  following  generalized  likeli- 
hood  function 


(1.3.2)  JI  /  Py  (r •  |a)  g(a;  8)  da 
i  G1  1 1  1 


should  be  maximized.  Since  this  method  coincides  with  the  maximum  likeli¬ 
hood  method  in  case  all  Ti  are  the  same,  we  called  it  the  general i zed 
maximum  likelihood  method.  A  preliminary  form  of  a  goodness  of  fit  test 


for  this  situation  is  discussed  in  Section  4. 


1.3.2.  BAYESIAN  ACCEPTANCE  PLANS.  Ir  this  section  we  assume  that 
the  prior  distribution  is  known.  We  further  assume  that  equipments  with 
values  of  the  parameter  a  in  some  subset  G0  of  G1  are  acceptable.  We 
confine  ourselves  to  acceptance  test  plans  of  the  form  (T*.  r*),  i.e.5  if 
in  test  time  T*  the  number  of  failures  of  the  equipment  is  less  than  or 
equal  to  r*,  the  equipment  is  accepted.  If  this  number  of  failures  is 
greater  than  r*,  the  equipment  is  rejected.  The  quantities 

p  -  Prob  {aeG0|r>r*} 

q  =  Prob  {a£GQ | r<r*} 


are  called  the  producer's  and  the  consumer's  risk  respectively.  These 
risks  are  usually  referred  to  as  posterior  risks  (cf.  Goel  and  Joglekar  (5)). 
These  risks  are  given  by  the  following  formul as 


(1  .3.3) 

(1.3.4) 


P  -  2 

r=r*+1 


r* 

q  •  z 

r=o 


Pj*  (r |a)  g(a;8)  da  j 
PT*(r jet)  g(a;8)  da  / 


2 

r=r*+l 


7Gi  pt*  (r  | a)  g(a;3) da 
Py*  (r|a)g(a;8)  da 


Once  the  consumer  decides  on  q  and  GQ,  he  can  determine  T*  and  r *  from 
(1.3.4)  (with  T*  smallest  possible)  and  thus  obtain  an  acceptance  test 
plan.  The  equation  (1.3.3)  will  then  determine  the  producer's  risk. 
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2.  EXPONENTIAL  TIME  TO  FAILURE 


2,1  THE  EXPONENTIAL  DISTRIBUTION.  The  most  widely  used  distribution 
for  the  time  to  failure  is  the  exponential  distribution.  This  distribution 
is  deemed  appropriate  for  electronic  equipment.  We  will  not  enter  into 
the  discussion  on  when  its  use  is  appropriate  and  when  its  use  constitutes 
a  misuse. 

The  density  function  is  given  by 

(2.1.1)  f(tje)  *  exp  (-t/0)  t>o,  9>o, 
the  reliability  function  is 

(2.1.2)  R( t | 0 )  a  exp  (-t/0) 

and  the  failure  rate  is  constant 

(2.1.3)  A(t)  =  1/9  t>o. 

It  follows  from  (1.1.1)  that  the  mean  time  to  failure  is  0,  i.e., 

(2.1.4)  9]  =  9. 

Because  of  the  strong  Markov  property  of  the  exponential  distribution 

(2.1.5)  Prob  {t>a+bjt>a}*  Prob  {t>b}  (a,  b>o), 

the  instantaneous  resurrection  and  the  instantaneous  renewal  assumptions 
lead  to  the  same  distributions  of  the  attribute  failure  data.  Thus, 
theorem  1.2.1  and  formula  (1.2.4)  yield 

(2.1.6)  PT(r | 9)  =  — Jq-  (T/0)r  exp  (-T/9), 
whereas  theorem  1.2.2  and  formula  (1.2.5)  yield 

(2.1.7)  f_{Tia)  ■  —  (Tr“ver)  exp  (-T/8) . 

r(r) 

Also  formula  (1.2.1)  is  reduced  to 

Prob  (5(t)  -  £(s)  =  r}  =  (^)r  exp  (- (t-s )/0 ) , 

i.e.,  the  stochastic  process  £(t)  is,  in  this  case,  a  homogeneous  Poisson 
process  with  constant  1/9.  The  formula  (1.2.3)  yields 

(2.1 .8)  9r  =  r9 

Because  of  (2.1.5)  the  mean  time  to  failure  equals  the  mean  time  between 
failures.  The  formula  (2.1.8)  tells  us  that  the  mean  time  to  the  rth 
failure  is  r  times  the  mean  time  to  failure. 
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2.2,  THE  TRUNCATED  INVERTED  GAMMA  PRIOR  DISTRIBUTION  AND  ITS  LIMIT  FORMS. 
The  parameter  of  the  exponential  distribution  is  the  mean  time  to  failure 
0.  In  the  set-up  of  section  1.3,  k=l ,  a=0,  G-.  *  (o,  +  »  ).  We  propose 
as  a  prior  distribution  the  truncated  inverted  gamma  distribution 

(2.2.1)  gO;Y,X,  a,  b)  -  YV  ilLil!  ^  exp  (-y/9)  .  3<ft<h 


where  Y>o  and 


r(x, Y/b)  -  r(x,Y/a) 


(2.2.2)  r  (X,  x)  ~  f  ux'Vu  du,  x>o. 


If  x=o  in  (2.2,2),  we  must  have  X>o  and,  in  this  case,  r(\,o)  =  r(X) 

If  X>o,  T(X,  x)  can  be  computed  either  by  numerical  integration  by  using 
the  formula 

(2.2.3)  T  (x,  x)  -r(x)  -  /x  uX-1  e"U  du 

o 

or  by  using  Pearson's  tables  (7),  where  the  function 

(2.2.4)  l(x,  p)  »  J_  /x  ^  UP  e'u  du 

r(p+i )  0 

is  tabulated.  From  (2.2.3)  and  (2.2.4)  we  obtain 

(2.2.5)  r(x,  X)  *  r(x)  (1-1  (x/A,  x-l)) 


The  recurrence  relation 


(2.2.6)  r(x+l ,  x)  *  xx  e'x  +  X  r(x,  x). 


the  asymptotic  formula 


e'"  (1+id  +  (M)  (*-2)  ,  .  .1 

x  x2 


as  x-*» 


as  well  as  the  limit 


o  if  x>l 
=  {  1/2  if  x=l 

1  if  -  . 27846... <x<l 


shed  some  light  into  the  behavior  of  this  function  (cf.  for  example 
Chapter  6  of  (8) ) . 

The  kth  moment  about  zero  of  the  distribution  (2.2.1)  is 


(2.2.7)  u/  =  yK  JAM 


r(X,  Y/b)  -  r(x,  Y/a) 


,  k=o,  1 ,  ... 
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We  now  look  at  some  interesting  limiting  forms  of  (2.2.1). 

Case  (A):  a-*o.  In  this  case  (2.2.1)  takes  the  form 

(2.2.1A)  g(9;  Y,  X,  o,  b)  -  yV<x+1>  exp  (-  y/9)  /V  (X,  Y/b) 

o<0<b,  y  >o 

and  its  moments  about  zero  are 

(2.2.7A)  uj*  =Yk  r  (X-k,  Y/b)  /  r(X,  Y/b),  k=o,  1,  ... 

Case  (B ) :  b-«°.  Here,  the  formulas  (2.2.1)  and  (2.2.7)  yield 

v  ^  )eXD  (  y/fti 

(2. 2. IB)  g ( 0 ;  y,  X,  a,  «  )  *  — - — * - - -  »  a<9,  Y  >o,  X>o 

T  (X)  -  r(x,  Y/a) 

(2.2.7B)  U'  »  y  k  IltJLL.  r(x"k»  Y/a) _ ,  0<  k  <X 

k  r(x)  -  r(x,  Y/a) 

Case  (AB):  a-*o  and  b-+®.  This  is  the  case  of  the  inverted  Gamma 
distri bution, 

(2. 2. TAB)  g < 9 ;  Y.X,  o,  •)  -  YV^X+1)  exp  (-Y/e)  /  r(X),  S>o,  y>o,  X>o 
(2.2.1AB)  ujj  -Yk/(X-1)  (X-2)  ...  (X-k),  o<k<X. 

Case  (r):  y+o.  In  this  case  we  obtain 

(2. 2. ID  g(0;  o,  X,  a,  b)  =  X9"(X+1)  /  (a*X-b'X),if  X^O  , 

(2.2.1 TA0 )  g(9;  o,  o,  a,  b)  =  1/9  log  (b/a),  o<a<9<b  . 

This  is  a  two-ways  truncated  Pareto  distribution.  Jf  X<o,  b  can  be  taken 
equal  to  +  •  in  (2. 2. IT),  whereas  if  X<o,  a  can  be  zero.  If  X=-l  (2. 2. IT) 
yields  the  uniform  distribution.  Thus 

(2.2.1FA_]_)  g(9;  o,-l,  a,  b)  =  l/(b-a),  o<a<0<b. 

If  b-*®  in  (2.2.irA_i)  we  obtain  an  improper  distribution. 

The  moments  are  easily  obtained: 

X.  (a-(A-k).b'('Vk))/(a’A-b'A)-if  k*X 
(2.2.7T)  u,-  -  {  x-k 

K  XX 

X  1 og ( b/a )  /  (a  -b  ), if  k=  X, 
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(2.2.7rAQ)  yj  =  ~~  (bk  -  ak)  /  log  (b/a),k=1,  2 . 

(2.2.7TA+B)  \i'k  =  Xak/(x  -  k),o<k<A, 

(2.2.7rA_A)  p'  =  bk  /  (k-A)  ,  k=o,  1 . 

(TA+B  means  that  in  2.2.7T,  Xis  positive  and  b-*»  whereas  TA_A 
means  that  Xis  neoative  and  a-+o.  The  uniform  case  is  included  in 
2.2.7?  with  X  *-l  j. 

2.2.1  ESTIMATION  OF  THE  PARAMETERS.  By  virtue  of  (2.1.6)  and  (2.2,1), 
the  marginal  distribution  (cf.  formula  (1.3.1)) has  the  form 

(2.2.8)  PT  (r;  Y,  X,  a,  b)  =  YX  Tr  [r(X+r,(T+f )/b)-r(x+r,(T+y)/a)]/ 

r!  (T+Y)X+r  [r(X,Y/b)  -  r(X,  v/b)] 


In  terms  of  the  incomple  Gamma  ratio  (the  function  given  by  (2.2.4)), 

(2.2.8)  takes  the  form 


PT(  v;  y,  x,  a,b)  =(T“)X  <^rY)  (A;r'"')  - 


.A+r-K  II 


)/b  vT7  ),+r- 


‘  1  |  I(Y/a/X,  X-l)  -  I  (y/b/X,  X-l) 

The  (descending)  factorial  moments  of  tnis  compound  distribution  are 

(2.2.9)  y(k)  a(T/Y)k  [r(X+k,y/b)-r(X+k,Y/a)]/[r(X,Y/b)  ~r(X,Y/a)] 

k=o,  1,  ... 

There  are  four  parameters  to  be  estimated.  If  we  have  a  sample  with  T 
fixed,  the  method  of  moments  can  be  used.  Let  m,. ,  denote  the  factorial 
sample  moments.  Thus,  by  matching  moments  we  get  ‘ 

(2.2.10)  m,^=(T/Y)k  [r(A+k;  Y/b)-r(X+k,Y/a)]/[r(X,Y/b)-r(X,Y/a] , lei  ,2,3,4 

By  dividing  the  kth  equation  in  (2.2.10)  by  the  (k-l)th  (remembering  that 
m(o)  3  u(o)  =  ^  anc*  us"in9  tbe  recursion  (2.2.6)  we  obtain 

(2.2.11)  m(k)/m(k_1)-(x+k-l)T/Y=(T/Y)[(Y/b);V+k'1e-y/b  -(Y/a)H'VY/a], 

[r(X+k-l,  Y/b)-?(X+k-l,  y/a)],  k=l ,2,3,4  . 

We  set  Mk  =  m^  -  (X+k-I)  T/Y  and  consecutively  divide  the 

the  equations  (2.2.11).  In  this  way,  we  obtain  the  equivalent  system 

fWrik  ^CWb/a)X+k  exDC(b-1-a'1)Y]/(l-(b/a)X+k_1  exp[(b'1-a”1  )Y1 ) 

k=l ,2,3, 
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(2.2.12) 

M1=(T/Y) (r/b )Pe"Y/b fl-(b/a)xexp  [(b^-a'^y] ]/ [r(A,y/b)-r(Wa) ] . 

A  little  more  algebra  with  the  first  three  of  (2.2.12)  yields 
abM3  -  (a+b)  TM£  +  T2^  =  0 

abM4  -  (a+b)  TM3  +  T2M2  =  C 

(2.2.13)  (b/a)A+1  exp  [Y(b'1-a’1)]  =  (bM2  -  TM- )/(aM2-TM^ ) 

The  first  two  of  these  yield 

ab  =  T2(M1M3  -  M2  )  /(H2  M4  -  M2  )=  T2B  (y,A,  T) 

a+b  =  T(M1M4  -  M2M3)/(M2M4  -  m|)  =  TA(y,A,T) 
which  Imply 

(2.2.14)  a  =  %T [A  -  (A2  -  4B)3*]  b  =  \  T  Ca+(aZ  -  AB)1/2]- 

Thus  we  have  a  and  b  exclusively  in  terms  of  v  and  A.  If  we  substitute 
these  values  of  a  and  b  in  (2.2.13)  and  the  fourth  of  (2.2.12),  we 
obtain  two  equations  in  y  and  A  which  can  be  solved  numerically. 

We  now  describe  a  simpler  estimation  method.  We  can  take  as  a  and  b 
the  lower  and  upper  confidence  bounds  for  9  as  estimated  from  the  failure 
data.  If  the  data  are  r . ,  i =1 ,  n,  and  if  r=L'r.  .then  e  =  nT/r  and 

I  rv  I 

a  two  sided  confidence  interval  for  e  at  confidence  level  1-ci  is 


(2.2.15) 


A 


2r9 


9 

(2r+2) 


A 

2  re 

x  a/2  ^ 


(cf.  p.  181  of  (2)). 

We  can  take  a  very  small  and  put  a  equal  to  the  lower  confidence  bound 
and  b  equal  to  the  upper  confidence  bound.  Then,  we  can  proceed  to 
the  computation  of  Y  and  A. 


It  is  easy  to  ..rite  down  the  maximum  likelihood  equations  or  the 
generalized  maximum  likelihood  equations  for  this  situation.  They  turn 
out  to  be  very  complicated  to  work  with  although  this  complication  is 
considerably  reduced  if  we  assume  that  a  and  b  are  known.  The  main 
disadvantage  is  that  numerical  integration  is  necessary  in  order  to 
compute  integrals  of  the  form 
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.X,  x;  =  /  UX-1  loo  u  e  u  du  -  r(Xrd»(X)-  / x  ux-l  logue-11  du. 
x  o 


However,  if  X  is  assumed  known,  this  difficuly  is  overcome.  If  X  and  y  are 
assumed  known,  the  situation  is  quite  manageable,  lie  omit  the  details. 

In  the  situation  of  type  II  censored  failure  data  we  obtain 

(2.2.16)  fr(T;  y,  X,  a,b)  =  y  py  (r»  Y»  X,  b) 

where  (r;  y,  X,  a,  b)  is  given  by  (2.2.8).  The  moments  about  zero  are 

(2.2.17)  u-  =  k  T(r+k)  r(X-K,  y/b)  -  r(X-k,  Y/a) 

k  Y  r(r)  r(x,  v/b)  -  r(x,  Y/a) 

An  analysis  similar  to  that  of  the  case  of  type  I  data  can  be  carried  out 

in  using  the  method  of  moments  to  estimate  the  parameters.  Incidentally, 

(2.2.16)  shows  that  the  generalized  likelihood  equations  will  have  the 
same  form  as  that  of  type  I  data  (since  r/j  does  not  depend  on  the 
parameters) . 

We  observe  that  the  posterior  distribution  is  again  a  truncated  inverted 
Gamma  distribution  with  parameters  T+y,  X+r,  a,  b,  i.e. 

hT  (9|r)  =  hr  (9|T)  =  (T+y)X+V(X+r+1 }  exp  [-(T+Y)/e]  * 


[r(X+r,  (T+y)/b)  -  T(X+r,  (T+y)/a)]. 

The  posterior  mean  is  given  by  (2.2.7)  for  k-1 ,  y  replaced  by  T+Y  and  X 
replaced  by  X+r. 

We  now  turn  to  the  more  hospitable  limiting  cases  of  the  truncated  inverted 
Gamma  prior  distribution. 

Case  (A):  a-+o.  The  marginal  distribution  for  type  I  attribute  data  is 
(2.2.8A)  PT(r;y,  X,  o,  b)=  YV  r(X+r,  (T+y/b)/r! (T+y)X+r  T(X,  Y/b) 

and  the  equating  of  sample  factorial  moments  with  the  theoretical  ones 
yields 

(2.2. 10A)  m(k)=  (T/Y)k  r(X+k,  Y/b)  /  T(X,  Y/b)  k=l,2,3. 

By  eliminating  b  we  obtain  the  quadratic  equation 

(2.2.18)  (ym^-U+UTm^))2  =  (ym,  ^ -XT)  (Ym^3p(,\+2)Tm^2^ ) 

from  which  y  can  be  expressed  as  a  simple  function  of  A, i.e.  y  -  y(x).  The 

parameter  b  is  given  by 
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(2.2.19)  b*T(Y  m(1)  -  ,\T)/(m(2)Y-  (A+l )  m(1) 


T)  5  b(A) 


Thus  we  have  both  y  and  b  exoressed  in  terms  of  X  and  we  only  have  to 
solve  the  equation  (first  of  (2.2.1 0 A ) ) 

(2.2.20)  [Y(X)fn(1}  -  XT]  r(X,  Y(x)/b(A)-T  [Y(x)/b(X)]X  exp(-Y(X)/b(X)=0 

This  can  be  done  by  usinq  interpolati ve  methods  and  Pearson's  tables  (7) 
for  the  values  of  r(X,  x)  (in  conjunction  with  formula  (2.2.5)).  We  omit 
the  details  for  type  II  data. 

Case  (B):  b-*».  In  this  case  the  method  of  moments  yields  the  following: 

The  quadratic  (2.2.18)  remains  unchanged,  in  (2.2.19)  b  is  replaced  by  a 
and  (2.2.20)  is  replaced  by 

(2.2.19)[Y(X)  m(])  -  XT]r(X,y(X)/a(X))+T  [y (X)/a(X)]Xexp(-Y(X)/a(X)=o  . 

Case  (AB):  a-*o  and  b-«°.  This  the  classical  case.  The  marginal  distribution 
for  type  I  data  is  a  negative  binomial  distribution.  The  moment  estimates 
as  well  as  the  maximum  likelihood  estimates  of  y  and  X  are  well  known.  In 
(6)  we  give  a  sufficient  condition  for  the  existence  of  the  generalized 
maximum  likelihood  estimates  of  y  and  X  when  the  test  time  for  each 
equipment  is  not  constant. 

Case  (T):  y+o.  The  marginal  distribution  for  type  I  data  (with 
repTacement)  is 

(2.2.81’)  PT(r;  o,  X,  a,b)  =  XT"X  [r(X+r,  T/b)  -  r(A+r,  T/a)]/r! (a"X-b"X) . 

where  A/(a~X-b*X)  must  be  replaced  by  1/log  (b/a)  if  A=o.  If  X+o,  -1,  -2,  -3, 
the  method  of  moments  yields 

(2.2.1  or)  m(k)  s —^iT  (a~U+k)-b"(X+l°)/(a'x-b~X),  k=l,2,3. 

From  these  equations  we  obtain 

a=  -J  [A-(A2  -  48)V2],  b=  [A  +  (A2  -  4B)  1/2] 

where 

A=[m(3)  X(X+3)  -  m(2)  (X+l)  (X+2 ) ]/  fm(1)m(3  )  U+1 )  U+3)-m22)  (A+2)2] 


S=[m^2^  A(X+2)-m^  (X+l}2]/[m^  (X+l)  (X+3)  -m^  ^+2)2]- 
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These  values  of  a  and  b  can  be  substituted  into  the  first  of  (2.2.1  Or ) 
(which  can  be  written  in  the  form 

b"(,)  -  »n/  [(x+i)  aM  .  xT] 

and  then  solve  numerically  for  X. 

If  X=o,  we  get  the  two  equations 

rn ^ -]  j  =  T(a-1  -  b_1)/log  (b/a),  tn^j  =  T2  (a"2  -  b"2 ) /2  log  (b/a) 
which  yield 

am(ij  (2m(2)  “  aTm(l })  l°g  (—^  -  1 )  =  2T  (m^  -  aTm^.) 

aTm(l ) 

b  =  2m^  —  3. 

If  A=1  (the  case  of  uniform  distribution),  we  get 

mO) 3  ~~  log  b/a>  m(2) 3  r2/ab 

which  can  easily  be  solved  numerically. 

If  X  *  -2  we  get 

m(1)  =  2T/(a+b) ,  m(2)  *  2T2  log  (b/a)/(b2  -  a2). 

Finally,  if  X  =  -3  we  get 

j  3  3T  (a+b)/2  (a2+ab+b2),  m^  =  3T2/(a2+ab+b2)  which 
imply  that  a  and  b  are  the  two  roots  of  the  quadratic  equation 

m(2)  ^  "  2m(l )  m(2)  T  X  +  (%(1)  ‘  ^(Z)5  jZ  =  °* 

Case  (TA):  y-*o,  a-o.  In  this  case  the  distribution  is  defined  if 
X<o  and  the  method  of  moments  is  applicable  only  if  X<-2.  Under 
this  assumption  we  have 

m ,  1  ^  =  XT/ (X+l )  b,  m^  *  XT2/(x-t°)  b~ 

15 


i 


which  imply  that 


A  =  1  '  C  m(2)/(m(2)  -  ]1/2*  b  3  ^T/(A+1 )  m(1) 


provided  that  m  <  m 


(1)  (2)  8  (1) 


Case  (TB):  y-+o,  b-*»,  X>o  .  In  this  case  we  get  the  estimates 


X  "  1+  Cm(2)/(m(2)  -  m^)i 


provided  that  m,  ,  >  nr 

(2)  (1) 


a  =  AT/  ( A+l ) 


Before  closing  this  section,  we  observe  that  the  question  of  identi fiabil ity 
of  Poisson  or  Gamma  mixtures  is  affirmatively  settled  by  methods  akin  to 
the  proof  of  uniqueness  of  Laplace  transforms. 

2.2.2.  ACCEPTANCE  PLANS  BASED  ON  POSTERIOR  RISKS.  Let  be  the 
minimum  acceptable  mean  time  to  failure.  In  the  notation  of  section  1.3.2, 
G0  ~  [9a  ,  b],  G^j  =  [a,  b]  and  the  equation  (1.3.4)  is  reduced  to 

^  ic  p  ic 

(2.2.22)  q  *  Z  PT*  (r;  Y ,  A,  a,  0O  )/  Z  PT*  (r;  y ,  A,  a,  b) 
r=o  r=o 

from  which  (r*,  T*)  with  r*  or  T*  minimum  can  be  computed  if  a  solution 
exists. 


2.2.3,  OTHER  PRIOR  DISTRIBUTIONS.  Numerous  types  of  distributions 
(e.g.  truncated  normal,  log-normal,  etc.)  can  be  tried  as  prior 
distributions  of  9.  There  is  an  extensive  literature  on  the  resulting 
Poisson  mixtures.  Because  of  lack  of  space,  we  will  not  deal  with  them 
in  this  report.  Most  of  the  information  can  be  found  in  the  encyclopaedic 
work  of  N.I.  Johnson  and  S.  Kotz  (9),  Chapter  3. 
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3.  THE  CASE  OF  NONCONSTANT  FAILURE  RATE. 


Under  the  instantaneous  resurrection  hypothesis,  which  is  usable  in  the  case 
of  repairable  aging  equipment,  the  whole  procedure  of  using  prior  distributions 
for  acceptance  plans  is  reduced  to  either  Poisson  mixtures  (theorem  1.2.1) 
or  to  transformed  Gamma  mixtures  (theorem  1.2.2)  just  as  in  the  case  of 
censored  tests  without  replacement  the  theory  of  Binomial  mixtures  (formula 
(1.2.6))  and  transformed  8eta  mixtures  (formula  (1.2.7))  is  the  relevant 
set-up. 

As  an  example,  consider  the  Weibull  time-to-fail ure  distribution 

(3.1)  f ( 1 1 9,  c)  *  c9_1  (t/Q)0"1  exp  (— (t/0 )c) ,  t>o,  0>o ,  oo. 

The  reliability  function  is 

(3.2)  R( t | 0 ,  c)  ^  exp  (-(t/e)c) 
and  the  failure  rate  is  given  by 

(3.3)  \(t)  =  C9"1  ( t/0  )c_1 . 


By  the  theorem  1.2.1,  the  distribution  of  the  type  I  censored  failure  data 
under  the  instantaneous  resurrection  assumption  is  given  by 

(3.4)  PT  (r J 0 ,  c)  =  (T/9)Cr  exp  (-(T/0)c) 

By  supposing  c  known  and  assuming  a  Gamma  prior  distribution  of  9"c, 
we  obtain  a  situation  almost  exactly  the  same  as  in  the  case  of  exponential 
time  to  failure.  The  difference  is  that  9  now  is  not  the  mean  time  to 
failure.  The  mean  time  to  failure  is  given  by 


(3.5)  91  =  9  r  (c'1  +  1). 


i.e.  there  is  a  change  by  a  factor  of  f(c-^  +1).  In  the  formula  for  the 
marginal  distribution  as  well  as  in  the  formulas  of  the  moment  estimators, 
maximum  likelihood  estimators  or  generalized  maximum  likelihood  estimators 


— _ 1 - - -*- 

v.naii^c  id  u  mc  i  cp  i  av-aucML  ui 


Mil  duaiuionai  minor  adjust- 


ment  has  to  be  made  to  the  formula  determining  the  acceptance  test  plans. 


Further  results  on  Bayesian  techniques  in  connection  with  the  Weibull 
distribution  can  be  found  in  (10). 
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4.  TOWARD  A  GOODNESS  OF  FIT  TEST  WHEN  THE  GENERALIZED  MAXIMUM  LIKELIHOOD 
ESTIMATORS  ARE  USED. 

When  the  available  attribute  failure  data  do  not  exhibit  uniform  time  of 
operation,  the  only  way  available  for  the  estimation  of  the  parameters  of 
the  prior  seems  to  be  the  generalized  maximum  likelihood  method  (cf.  section 
1.3.1,  formula  (1.3.2)).  In  order  to  test  goodness  of  fit  in  this  case, 
the  ordinary  y}  will  not  do  simply  because  the  usual  statistic  is  not 
asymptotically  y}  -  distributed  in  this  case. 

Suppose  that  we  have  decided  the  ranges  of  the  cells  and  let  C-j ,  Co,  ....  Cj, 
be  these  cells  (e.g.,  suppose  that  all  equipments,  which  had  no  failures 
belong  to  cell  Ci,  the  equipments  which  had  one  or  two  failures  belong  to 
cell  C2  etc.).  Suppose  further  that  the  sample  size  is  n  and  that  the  ith 
equipment  had  r^  failures  in  time  T^ ,  i*l,  2,  ...  n.  Let  P.^  be  the 
probability  that  the  ith  equipment  belongs  to  the  cell  Cj,  i.e. 

Pij  3  reCj  Pt1  *r)  ls1,  2 .  n  5  js1,  k 

Obviously, 

k 

(4.2)  I  P,,  =  1. 

3*1  J 

We  set 

(4-3)  pj  1 4  j,  V 

Then,  obviously 

(4.4)  Z  Pj  a  1 , 

j=l 

and  the  expected  cell  frequencies  are  nPj,  j=l,  ....  k. 

Let  yu  j*l,  k  be  the  observed  cell  frequencies.  The  joint  moment 
generating  function  of  y-j ,  . ..,  yk  is 

n  k  t . 

H  Z  Ptje  J. 
i=l  j-1  J 

Therefore,  the  joint  moment  generating  function  of  the  random  variables 

Xj  =  yj  '  15  k  n  k  . 

(4.5)  \p  (t],  ....  tk)  =  exp  (-n  Z  Pi  t.)  n  Z  P-jje1^ 

j-1  J  J  1=1  j=l 

From  the  McLaurin  expension  of  this  function  we  obtain  the  following  joint 
moments  that  are  of  interest  to  us  (the  indices  range  from  1  to  k). 
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(4-6)  maa  =  nPa  -  Z  .ja,  mab  =  -  Z  p i a p i b »  (afb) 


(4*7)  maaaa  =  3n^a  +  npa  '  6nPa  zP?a  ’  7zpia  +  12  2pia  + 

2  2  4 

+3{ZPia)  -  6ZPia  2 

(4.8)  maabb  ”  n^papb  “  ZPiapib  "  nPa  2Pib  "  nPb  2Pia  +  2  ZPiaPib  + 

*  2Epfapjb  +  <»?«>  <2Pfb)  +  2<ZPiaPib>2  ’  6EP1aP!b  ' 


We  consider  the  statistic 

2  n 
(4.8)  S  =  Z  x  ^/wi ,  Wa  =  (nP.-  -  Z 
J  J  J  J  i=1 


2 

PijO/O-Pj) 


The  mean  of  this  statistic  is  given  by 

k  k  k 

(4.9)  E(S)  =  Z  m,,/w,  =  Z  (1-P,-)  =  k  -  Z 


j-1 


k-1 


j-1 


j=l 


where  (4.6),  (4.8)  and  (4.4)  have  been  used.  The  variance  is  given  by 
(4.10)  Var(S)  3  E(S2)  -  E(S)2  =  Z  ~  2 


.  maaaa/w, 
a=1  a 


+  2  Z 


a^b  rt'aabb/wawb  ^ 


where  the  values  of  m 
The  Chebyshev  inequal 
statement  about  the  goodness  of  fit. 


a  and  maabb  are  given  by  (4.7)  and  (4.8)  respectively, 
ity  can  be  used  to  obtain  a  (weak)  probabilistic 


If  we  adjust  4.5  because  of  the  weights  w.-  (i.e.,  if  we  consider  the  joint 

moment  generating  function  of  the  variables  ,  jal  4 _ .  k)  take  its 

logarithm  and  expand  into  Taylor  series  aboutJtheJorigin,  we  obtain  that 
for  n-*-  ”,  the  variables  are  asympotically  normally  distributed  with 
variance-covariance  matrix  of  rank  k-1 


The  non-zero  eigenvalues  of  this  matrix  are  not  all  equal  to  1  (in  general) 
and  therefore  the  asymptotic  distribution  of  the  statistic  S  is  not  a  x2 
distribution.  The  cumulants  of  the  asymptotic  distribution  of  S  are  given 
by  (cf.  Chapter  29,  Section  3  of  (11)). 

(4.11)  «cm  =  2m_1  (m-1 )  !  t r ( Am )  m=l ,  2,  ... 

Thus,  in  particular 

=  E(S)  =  tr(A)  =  k-1 
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2 

<2  =  w2  =  Var( S)  =  2tr  (A  ) 

<3  *  u3  3  8tr(A^) 
k  »  -  3^2  =  48tr(A^)  etc. 

It  is  sometimes  more  convenient  to  approximate  the  variance  of  S  this  way 
than  compute  it  through  (4.10). 

We  did  not  study  the  effect  on  S  of  the  estimating  of  parameters  from  the 
data  used  to  form  S.  Hopefully  something  not  very  far  from  the  classical 
case  of  "reduction  of  degrees  of  freedom"  will  occur. 
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